AN EXISTENCE RESULT FOR THE SANDPILE PROBLEM 
ON FLAT TABLES WITH WALLS 



GRAZIANO CRASTA AND STEFANO FINZI VITA 

Abstract. We derive an existence result for solutions of a differential system which 
characterizes the equilibria of a particular model in granular matter theory, the so-called 
partially open table problem for growing sandpiles. Such result generalizes a recent 
theorem of [6] established for the totally open table problem. Here, due to the presence 
of walls at the boundary, the surface flow density at the equilibrium may result no more 
continuous nor bounded, and its explicit mathematical characterization is obtained by 
domain decomposition techniques. At the same time we show how these solutions can be 
numerically computed as stationary solutions of a dynamical two-layer model for growing 
sandpiles and we present the results of some simulations. 



1. Introduction 

In the last years an increasing attention has been devoted towards the study of differen- 
tial models in granular matter theory (see, e.g., [2] for an overview of different theoretical 
approaches and models). This field of research, which is of course of strong relevance in 
the applications, has also been the source of many new and challenging problems in the 
theory of partial differential equations (see, e.g., [3, 6, 11, 16]). 

In this paper we deal with the rather simple phenomenon of the evolution of a sandpile 
created by pouring dry matter on a fiat bounded table. In such a model the table is 
represented by a bounded domain C M^, and the time-independent vertical matter 
source by a nonnegative function / G L^(J7). Recently, Hadeler and Kuttler [15] proposed 
a new model, extending the ones studied in [4] and [5], where the description of the heap 
evolution is based on the observation that granular matter forms heaps and slopes (the 
so-called standing layer), while small amounts of matter move down along the slopes, 
forming the so-called rolling layer. We also mention [16], where Prigozhin has studied, 
both from the theoretical and numerical points of view, a degenerate parabolic problem 
and its equivalent formulation as a variational inequality, and [3] , where a similar approach 
has been used for growing sandpiles on the whole plane. It is worth to remark that the two 
different dynamical models of [15] and [16] (see for example [17] for a comparison between 
them) have theoretically the same set of admissible equilibria. 

Let us denote by u(t, x) and v{t, x), t >0, x G fl, respectively the heights of the standing 
and rolling layers. Neglecting wind effects, the local dynamics depends on the local slope 
\Du\ and on the local density of rolling matter v (that we shall call transport density). For 
stability reasons, at any equilibrium configuration the local slope \Du\ cannot exceed a 
fixed constant (the critical slope), that we normalize to the value 1, and it must be maximal 
where transport occurs (that is, where v > 0). Moreover, we assume that the matter falls 
down the table when the base of the heap touches a portion T of the boundary of f2. Prom 
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the model point of view we are thus assuming that on dii \ F we have a (arbitrarily high) 
vertical wall, while on F the table is "open" . In the following, we shall refer to the open 
table problem in the case T = dil, whereas in the partially open table problem, T will be 
a nonempty closed subset of dO,. 

The dynamical model proposed in [15] deals with the open table problem, but it can be 
extended to the partially open table problem in the following way: 

dtv = dw{v Du) - (1 - \Du\)v + / in [0, oo) x O 
dfU = (1 - IZ^^D^; in [0,oo) x Q 

(1) { u(0, •) = v{0,-) = inn 

du 

u = on r , V— = on dn\T , 
ou 

where ^ denotes the normal derivative of u. The nonlinear term which appears in the 
previous equations with opposite signs represents the exchange term between the two layers 
during the growth process. Before the equilibrium has been reached, a partial surface flow 
is allowed also at sub-critical slopes. As far as we know, a rigorous theory for this model is 
not known, and its equilibrium configurations have not been characterized in the general 
case. For the open table case, a finite difference scheme has been proposed in [13], which 
offers at the same time a tool for the numerical description of stationary solutions. 

From the physical considerations above, an equilibrium configuration (n, v) for (1), with 
u,v nonnegative functions in must satisfy 

— div(i; Du) = f infl, 
\Du\ < 1 a.e. in 



(2) 



\Du\ = 1 in {i; > 0}, 

= on F, V — =Oondn\T. 

ou 



In the case of the open table problem (i.e., F = solutions of (2) have been completely 
characterized by Cannarsa and Cardaliaguet [6] (see also [7] for an extension to higher 
dimensions). More precisely, denoting by dg^: O — > M the distance function from the 
boundary of fi, they proved that there exists a nonnegative continuous function Vf (with 
an explicit integral representation) such that {dQ^,Vf) is a solution of (2), and any other 
solution {u,v) satisfies v = Vf in Vt and u = dgQ in {vf > 0}. 

Aim of this paper is to extend these results to the case of the partially open table 
problem. As we shall see in the sequel, the presence of vertical walls has a relevant influence 
on the regularity of stationary solutions. Namely, we cannot expect the transport density 
i; to be a continuous function as in the case of the open table problem. This fact has 
several consequences. 

First of all, we cannot give a pointwise meaning to the boundary conditions in (2). Our 
choice here it to set v in L^{Q,), and to consider a weak formulation of the problem (see 
(9) below). Another possible choice, which we do not pursue here, could be to set v in 
the class of functions with bounded variation, so that the trace of v on dfl is well defined. 
Our main result (see Theorem 2.2 below) states that there exists a nonnegative function 
Vf € L^(Jl), of which we give an explicit integral representation, such that {dr,Vf) is a 
solution to the weak formulation of the problem, where dp denotes the distance function 
from F (see (5) for its precise definition). 

The second major consequence of the lack of continuity of v concerns the uniqueness of 
solutions. Namely, the uniqueness of the transport density v for the open table problem 
was proved in [6] using a blow-up argument, first introduced in [11] for the analysis of 
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mass transport problems in the framework of the Monge-Kantorovich theory, that rehes 
on the continuity of v. In our opinion this argument cannot be adapted to the case v e L^. 
Nevertheless, it could be possible to prove a uniqueness result in a restricted class of more 
regular functions c L^{^), provided that one can show that Vf e J^. 

A precise formulation of the existence result mentioned above requires some notation. 
The table O C will be a Lipschitz domain, i.e. an open bounded connected set with 
Lipschitz boundary, and the open boundary F C dQ will be a nonempty closed subset of 
d^}. Let d: CI X CI ^ [0, oo) denote the path distance in CI, defined by 

(3) d{x,y) = inf{length(7); 7: [0, 1] — O Lipschitz path joining x to y} , 
and let 

(4) Lip^(n) = {u: f2 — M; u Lipschitz, ^(a;) — u{y) < d{x,y) Vx,y G CI} 

be the set of 1-Lipschitz functions in CI with respect to the path-metric d. Let us denote 
by dr' CI ^ [0, 00) the path distance function from F, defined by 

(5) dr{x) = inf d{x,y), x eQ. 

It is easily seen that, if F = dCl, then dr is the Euclidean distance from the boundary of 
CI. Moreover, dr belongs to the space of functions 

(6) Lipf (O) = {ue Lip^O); M = on F} . 

It is also known that dr is the maximal function among all functions in Lipp(r2). Thus dr 
is the maximal size of the standing layer. 

Given a nonnegative function / G L^{Cl), let uf. Cl —> [0, cxd) denote the function defined 

by 

(7) '^f{x) = m.ax{Gz{x); z G supp(/)} x & Cl , 
where, for every x,z G Cl, 

dr{z) — d{z, x), if d{z, x) < dr{z), 
0, otherwise , 



(8) Gz{x) 



and supp(/) denotes the (essential) support of /, that is, the complement in Cl of the 
union of all relatively open subsets A CCl such that / = a.e. in A. 

It is readily seen that the graph of the positive part of is, in the path metric, the 
maximal cone of unitary slope, with apex in z, whose base is contained in Cl and touches 
F at some point. It is clear that < Uf < dr, and that Uf £ Lip^(r2). Moreover, since Uf 
is the sup-envelope of all the cones Gz with z G supp(/), it is plain that Uf represents the 
minimal standing layer for an equilibrium configuration with respect to the given support 
of the source. 

We can now state the weak formulation of problem (2): Find nonnegative functions 
u,v: Cl ^ [0, +00) such that 

{u G Lipp(i7), V G L^iCl), u,v >0, 

I vDu .D<P=[f<P V(/) G C~(M2 \ r) , 
Jn Jn 

\Du\ = 1 a.e. in {v > 0}, 

where C^(M? \ F) denotes the space of C°° functions R with compact support 

in \ F. In this weak formulation the boundary conditions on u and v are embedded in 
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the choice of the test functions space C^{M? \ T) and in the condition u G Lipp(ri). It is 
readily seen that, if u and v are smooth enough, tlicn problem (9) is equivalent to (2). 

The plan of the paper is the following. In Section 2 we state the main existence result, 
constructing explicitly a transport density Vf (see formula (16) below) and showing that 
a pair {u,Vf) is a solution to problem (9) for every u G Lipf(f2) satisfying Uf < u < dr- 
Moreover, we show that no other function u G Lipp(J7) can be the standing layer for the 
problem. Section 3 contains the proofs of these results, that are mainly based on a Change 
of Variables formula of some independent interest (see Theorem 3.3 below). In Section 4 
we compute the explicit solution in a simple case, which will be compared in Section 5 to 
the numerical equilibrium solution of the dynamical model (1) obtained via some finite 
difference schemes. Indeed, in this last section the specific difficulties of such a numerical 
approximation are discussed in details. 

2. Existence of a solution 

Throughout this paper, C and T C d^} satisfy the following assumptions. 

(HI) O C is a Lipschitz domain, i.e. a nonempty open bounded connected set with 
Lipschitz boundary. 

(H2) r = [J^iTi is a nonempty closed subset of d^l, with Fi, . . . ,rAr connected arcs 
of 9ri, pairwise disjoint (up to the endpoints) and of class C^. We denote by Ai, 
Bi the endpoints of the arc Fj, z = 1, . . . , A'", and by the collection of all these 
endpoints. 

(H3) For every a; G there exists y G F such that dY{x) = \x — y\. 

Remark 2.1. Observe that (H2) does not prevent the intersection of two arcs at the end- 
points. This is the case, for example, if J7 is a square and Fj, i = 1,...,4, its sides. 
Nevertheless, in order to simplify the proofs, in the following we assume that the arcs Fj 
do not intersect at the endpoints. The general case can be treated with minor modifica- 
tions. 

Condition (H3) says that, for every a; G fi, there is a point y G F such that the closed 

segment [x, y] with endpoints x and y is a path of minimal length joining x to F. 
There are three relevant cases where these conditions are satisfied: 

(a) F = do., and is a domain of class C^. 

(b) F = 50, and O is a domain with piecewise boundary. 

(c) n is a non-empty open bounded convex set, and F satisfies condition (H2) above. 

Case (a) and (b) refer to the open table problem. They were considered respectively in 
[6] (see also [7] for an extension to ri C M") and [14] in the case of piecewise C^'^ boundary 
with outer (i.e. "convex") corners. 

Let F* = F \ F^. For every y G F* let v{y) denote the inward unit normal vector of 50, 
and let n{y) denote the curvature of at y. 

For every x G we denote by 

nr(x) = {y G F; dr{x) = d{x,y)} 
the set of all projections of x on F. By (H3) it is clear that 

nr(a;) = {y G F; dr{x) = \x- y|} . 
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Figure 1. Open boundary T (dotted), extended ridge TZ (bold) and trans- 
port rays 

For every x € and every y £ nr(x), let us define 

(10) l{x) = dr{x) ■ sup{t > 0; y + t{x - y) £ Q and y € Uriy + t{x - y))}, 

x-y 



(11) 
(12) 



m(x) = y + Ux) 

\x - y\ 

Rx = }y,m{x)l, 



where \x, y\ denotes the segment joining x to y without the endpoints. The set Rx is called 
a distance ray (or transport ray) through x and, in general, it depends on the projection 
point y. On the other hand, it is not difficult to prove that m(x) and l{x) do not depend 
on the choice of y G nr(a;). It is readily seen that l{x) is the length of the transport ray 
Rx- The set 

(13) n = {x£U; drix) = l{x)} 

will be called the extended ridge of O (see Figure 1). It coincides with the usual definition 
of ridge when F = (see [12]). Finally let r: $7 — > [0,oo) denote the normal distance to 
the ridge, defined by 

(14) t{x) = l{x) - dr{x), xgU. 

This function is clearly bounded from above by max{dY{x); x € Q}. 

Let $7* C be the set of regular points of dr , that is the set of all x £ 0, such that the 
projection nr(a;) is a singleton. It is well-known that £^(^7 \ 17*) = 0. 

Let x G , let nr(x) = {y}, and define for every t £ [0, t(x)] 

drix)+t 



(15) 



Mxit) 



drix) 
l-{dr{x)+t)K{y) 

1 - dr{x)K{y) 



if y G F^ 
if y G F* . 



Let us define the function 



fT{x) 



(16) 



Vf{x) 



f{x + tDdr{x))Mx{t)dt, ifxGJ7*, 
0, if X G S7 \ 17* . 

We remark that, if ]y, z[ is a transport ray, then z £ TZ, and 

lim Vfix) = . 

x^z, x^}y,zl 

In other words, the transport density Vf vanishes at the end of each transport ray. 
The main theoretical contribution of this paper is the following existence result. 
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Figure 2. An example of domain decomposition 



Theorem 2.2 (Existence). Let 0, and T satisfy (H1)-(H3). Then the pair (u,Vf) is a 
solution of (9) for every u G Lip\{^) satisfying Uf < u < dr- Furthermore, if {u,v) is 
any solution of (9), then Uf <u< dr- 

It is apparent that the result above characterizes all the possible standing layers as the 
functions u E Lipp(J7) satisfying Uf < u < dr- On the other hand, the uniqueness of the 
transport density vj remains an open problem. 



3. Proof of Theorem 2.2 

In what follows we shall always use the following decomposition of the set 0,* of regular 
points of dr- For every i = 1, . . . , A^, let us define the sets 

n* = {xen*; Ur{x) er*}, 

(17) 

nf = {xe n*; Ur{x) = {Ai}}, nf = {x € n*; Ur{x) = {Bi}} . 



Then (see Fig. 2 for an example) 



N 

B 



1=1 

Next, we need an approximation of this decomposition in terms of sets that can be 
easily parameterized (see Theorem 3.2 below). For every e > let us define the sets 

(18) G^ = (xGM^; min|x-y| <el, Tj = dGj, i = l,---,N. 

From Remark 2.1 and (113), there exists r > with the following property: for every e G 
(0, r], the sets Gj are pairwise disjoint and of class C^'^- Given e E (0, r] and i = 1, - - - , N, 
for every y G F- let I'^^y) denote the outer (with respect to Gl) unit normal vector to F- 
at y- In the following, e will always denote a number in (0,r]. 
Let us define the maps 

(19) ci>^^: F^ xM^m2, <^l(y,t) = y + {t-e)u'{y), y G F,S t G M, i = l,...,iV. 

The following lemma states that 0, can be parameterized using the maps ^j- Moreover, 
this parametrization is independent of the choice of e. Let us define 

Ti = Ttnn, i = i,.--,N, F^ = Un, r^ = Un, 

i=l 1=1 

and denote by 11*^ = Ur^ the projection operator on F"^. 
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Figure 3. The set 



Lemma 3.1. Let x G i7 and let y E nr(x). If y & Ti, then for every e S (0, r) there exists 
a unique ye G such that y^ S W{x) and x = ^l{ye,dr{x)). Moreover, x G ]y,ye[; if 
drix) < e, and G G f,^ if drix) > e. 

Proof. Let y^ be the intersection of L^ with the ray starting from y and passing through 
X. From the definition of F^ it is straightforward to check that y^ satisfies ah the stated 
properties. □ 

Let j/g € F*^, and let y G F be the unique projection of on F. Let us define 

'-1/e, if^/GF^ 
K{y)/{l-eK{y)), if y G F* . 

It can be checked that K'^(ye) is the curvature of F*^ at every point y^ where the curvature 
is defined. We remark that, for every z = 1, . . . , A^, k'^ is defined at all but four points of 
f^. Let X G O*, let W{x) = {yj, and define 

l-{dr{x) + t-e)K'{y,) 



(20) 



(21) 



Mm 



t G [0,r(x)]. 



1 - (dr(x) - e)K^(y,) ' 

It is easy to check that, if nr(x) = {y}, then 

(22) Ml{t) = MS). VtG[0,T(x)]. 

The main tool needed for the proof of Theorem 2.2 is a Change of Variables formula 
(see Theorem 3.3 below). As a first step, we need the following approximation result (see 
Figure 3). 

Theorem 3.2. For every e G (0, r] and i = 1, . . . , N , let 



N 
i=l 



QH{$,^(y,t); yGF^ 0<t</(y)}, 

Then C n, and lim.^o \ = 0. 
Proof. It is not difficult to see that the set 

N 

0':= \J{^t{y,ty, yGF^ 0<t<liy)} 

i=l 

covers f2 up to a set of Lebesgue measure A^vre^, i.e. £^(0 \ O^) < Nire^. Moreover, from 
Lemma 3.1 it is clear that C C The conclusion now follows observing that 

N 

0^\n^QD:= U{cI>Ky,/(y)); yG^}on, 
1=1 
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and that D has vanishing Lebesgue measure (see [9], Corollary 6.8). 

Let us define the map : x M — by setting t) = t) ii y e f|. 
Theorem 3.3 (Change of Variables). For every h G L^{Cl) we have 



(23) 



r r \ rKv) 

/ h{x) dx = lim / / h{^'{y, t)) [I - {t - e)K,'{y)\ dt 
Jn <^-^o Jr^ Jo 



lim 



h{y + tv'{y))[\-tK'{y)\dt 



dn\y) 
dn\y). 



Proof. Prom Theorem 3.2 we have that 



N 

h(x) dx = lim / h(x) dx = lim / h(x) dx . 

1=1 I 



For every z = 1 , . . . , we have the decomposition 

= nr u n^^ u nf , 

defined as in (17) replacing il* with $7*^. 

Using the arguments developed in [9, Sect. 7], it can be checked that 



(24) f h{x)dx= L [ ' h{^^{y,t))[l-{t-eX{y)]^ 
We shall now prove that 

(25) / ^ ^ h{x) dx= i ^ ^ f^'^ h{^'{y, t)) [1 - (t - e)Hf{y)] 



dt 



dn\y). 



dn\y). 



Let us define the region 



A, = {^'{y, t)■y(^^f^ nl'"", 0<t< l{y)} 



It is clear that C ^l'^- Hence, it is enough to prove that 



(26) 



h{x) dx 



liy) 



h{^'{y,t))[l-{t-€X{y)] dt 



dn\y) 



□ 



Observe that, by the very definition of k^, we have that nf'iy) = — 1/e for every y G V\r\Q.f. 
Then, the integral in brackets becomes 



/ hi^'iy,t))[l-it-e)^'iy)]dt= h{^'{y,t)) - dt . 
Jo Jo e 

On the other hand, the integral over A^ can be computed using polar coordinates {p,9). 

We remark that n fll'^ is an arc of circumference of radius e, so that dTi}{y) = edO, 

hence formula (26) follows. 

It is clear that (25) holds if O^'"^ is replaced by ^l'^ ■ The identity (23) then follows 

from (24) and (25). □ 

Lemma 3.4. The function Uf defined in (7) satisfies Uf = dr in the set {vf > 0}. 
Moreover, if u & Lip^{il) is a function satisfying Uf < u < dr in il, then Du = Dd^ in 
{vf > 0}. 



Kv) 
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Proof. From the very definition of Uf it is plain that Uf{x) > dr{x) for every x G supp(/). 
On the other hand, by the maximahty of c?r among all functions of Lip^(ri) vanishing on 
r, we conclude that 

(27) Uf = dT in supp(/). 

Now let X E {vf > 0} and let us prove that Uf{x) = dr(x). From the definition of 
we have that x G 17*, that is, nr(.T) is a singleton {y}, with y G F. Let [y, -z] be the 
transport ray through x. Since f/(x) > 0, we have that Jx, 2;|n supp(/) / 0. Let xq be 
a point belonging to this intersection. Prom (27) we have that Uf{xo) = dr{xo), which in 
turn implies that Uf = dr along the segment [y,.xo]. Since x G |y,a;o]], we conclude that 
Uf{x) = dr{x), and the first part of the lemma is proved. 

The second part follows from [1, Prop. 4.2] (see also [8, Lemma 7.3] and [10, Lemma 4.3]) 
upon observing that u = dr in {vf > 0} . □ 

Lemma 3.5. Let {u,v) be a solution of (9). Then 

(28) j /ii = max|y" fw; w E Lipl{Vt)^ . 
Proof. By a standard approximation argument we have that 

(29) / vDu-D4>= [ f(p V0: n^M Lipschitz, = on P. 
Ja Jn 

Let w G Lipf (n). Since v >0, \Dw\ < 1, \Du\ < 1, and \Du\ = 1 a.e. in {v > 0}, we have 
that 

[ vDu- {Dw - Du) < . 
Ja 

Hence we infer that 

[ fu - I fw> I vDu- {Dw - Du) - f{w -u)=0, 
Jn Jn Jn 

where the last equality follows from (29) with (j) = w — u. Since this inequality holds for 

every w G Lipp(i7), (28) follows. □ 

Proof of Theorem 2.2. We divide the proof into three steps. In the first one we shall prove 
that {dr,Vf) is a solution to (9). Then, in the second step, we shall prove that also {u, Vf) 
is a solution for every u G Lip^(r2) satisfying Uf < u < dr- Finally, in the last step we 
shall prove that if m G Lipp(O) is a nonnegative function with {u < Uf} ^ 0, then for 
every nonnegative function v G Li^ip.) the pair {u^v) is not a solution of (9). 

Step 1. We give only a sketch of the proof, since it follows the lines of the proof of 
Theorem 7.2 in [9]. 

Given 4> £ C^(M^ \F), we have to prove that (9) holds with u = dr- Using the Change 
of Variables formula (23) we have that 

(30) //</. = lim/ (^''\{^\y,t))f{^\y,t))\\-{t-e)K\y)\dt dn\y) . 

For every fixed y G P*^, let us integrate by parts the term in brackets. Recalling that 
$•^(^,0) = y — ev'^iy) G P, we have that 4'{^'^{y,0)) = 0, hence the integration by parts 
gives 



I{y) := j " 4>{^\y,t))f{^\y,t))[l-{t-e)K^{y)\dt 
Jo 

= f^'^ D<f>{^'{y, t)) ■ v'{y) f^"^ /($^(y, s)) [1 - (s - e)K^(y)] ds dt . 
Jo Jt 
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From the definition (16) of Vf and (22) we deduce that 

- (i - e)K^{y) 



ds , 



so that 



I{y) = f^"^ D<P{^'{y,t)) ■ v\y) Vf{<^\y,t)) [I - {t - e)K\y)\ 
Jo 



dt. 



Observe now that DdY{^^{y,t)) = v^{y) for every y G f ^ and every t G (0, Z(y)). Substi- 
tuting the last expression for I{y) in (30) and using again the Change of Variables formula 
(23) we finally get 

/ fcl> = \im L f^''\vfDcl>-Ddr){^%y,t))[l-{t-eX{y)]dt dH\y) 

= J VfD(j)-DdT. 

By the way, choosing (f) = dr (see (29)) we have that 

/ vf = fdr< +00, 
Jn Jci 

hence the nonnegative function Vf belongs to L^{Q). 

Step 2. Let u G Lip-'^(n) satisfy Uf < u < dr- Prom Lemma 3.4 we have that u = dr and 
Du = Ddr in the set {vf > 0}, hence 



VfDu-D(j)= I VfDdvD(i)= I f4> \/(t>€C^ 



\r), 



where the last equality follows from Step 1. 

Step 3. Let u G Lipp(ri) be a nonnegative function with {u < Uf} / 0. From [10], 
Proposition 4.4(iii), we deduce that also the set {x G supp(/); u{x) < dr{x)} is not 
empty and, in particular, u < dr on sl set of positive measure where / > 0. Since u < dr 
and / > 0, we infer that 



(31) 



[ fu< [ fdr. 
Jn Jn 



On the other hand {u,v) is a solution to (9), so that (28) holds, in contradiction with 
(31). □ 



4. A TEST EXAMPLE 

In this section we describe a simple example which illustrates very well how the presence 
of vertical walls on the boundary can influence the regularity of solutions (u, v) of (9). Let 
n = (0, if be the unit square of M^, and F = {{x, y) G : < x < 0.5 ; ?/ = 0} the only 
open part of its boundary. Assume / = 1 in all of O. Prom the picture in Pig. 4 we see 
that the sand transport rays behave differently in the two half sides of the table: in the 
left-hand side they lay parallel in the direction of F, whereas in the right-hand side they 
converge all together into the extremal point P = (0.5, 0), creating a singularity. 

Since / = 1 in fi, we have that Uf = dr so that the only possible standing layer is 
u = dr- The explicit computation for the solution {dr,Vf) of Theorem 2.2 can be done 
by decomposition of the domain 0, along the segment PQ. Using polar coordinates (r, 6) 
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Figure 4. The domain 0, of the example: P is a singular boundary point, 
the solution v results discontinuous along the line PQ 



centered in P (with 6 £ [0, |]) in the right hand side, from (15) and (16) we get in 
particular 

( 1-y, if X < 0.5, 

(32) Vf{x,y) = \ 

[Jl^'^^.dp, ifx>0.5, 

where l{6) denotes the length of the transport ray from P to the ridge on the wall boundary 
along the 9 direction (see (10)). It results that f/ G L^{0.) is unbounded near P, it is 
discontinuous along the segment PQ, and its gradient is discontinuous along the segment 
PS. The graph of the functions dr, Vf and their level lines are shown in Fig. 5. 

5. Numerical detection of stationary solutions 

In [13] the numerical approximation of the two-layer model of [15] was studied to sim- 
ulate growing sandpiles on an open flat table. Here we have considered the natural gener- 
alization (1) of such a model in the case of the partially open table problem, in order to 
get solutions of (9) as equilibrium solutions of a system of two evolutive partial differential 
equations. The extension of the finite difference scheme introduced in [13] to such a system 
is enough straightforward. For a given discretization step h = Ax, we introduce in the 
domain 0, (for simplicity, a rectangle) a uniform grid of nodes Xij, and we denote as usual 
by (tifj,v"j) the components of the discrete solutions at time = nAt. Then our fully 
explicit finite difference scheme can be written as 



(33) 
(34) 
(35) 
(36) 
(37) 



n 



,n+l 



D<3 



+ At(l 



yi,j, 

Vn, 

if Xij € O \ r, Vn, 



if Xij £ r 







where the discrete gradient vectors Du" and Du" are computed respectively, component 
by component, through the maxmod and the upwind finite difference operators, and D^u 
denotes the standard five-points discretization of the Laplace operator on the grid (see 
[13] for the details). What is new in this scheme is the wall boundary condition (37), 
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20 40 60 80 100 20 40 60 80 100 

X X 



Figure 5. Exact solution {dr,Vf) of the test example with level lines 

whose implementation requires some comments. The standard way is the following: after 
(33) and (34) have been applied, we look for the sign of (Du"'~^^ ■ v) at the wall nodes. 
If it is strictly positive (as it happens for nodes which are in the extended ridge TZ, that 
is which are starting points of a transport ray to F) then v^^^ is set to zero. If this is 
not the case, one should modify u^~j^ on the boundary in order to fulfill (37). This is 
not the best strategy. In fact this situation corresponds to the pathological case of nodes 
belonging to boundary transport rays, that is when there exist straight portions of the wall 
boundary, as in the test example of Section 4. Referring to Fig. 4, on the west side of the 
square the sand flow is parallel to the boundary (and also to the mesh in this particular 
case) and there is no need to impose any boundary condition: the discrete solution 
naturally satisfies a no flux condition at those points. Also the south portion of the wall 
in the example coincides with a transport ray, but in that case the normal derivative of 
u is naturally negative in the boundary nodes, and it becomes zero only asymptotically 
in time (at the equilibrium). Then, by continuity arguments, the best choice seems to us 
simply to impose a no flux boundary condition for v at those points. 

The direct application of scheme (33)-(37) is anyway not so efficient, due to the numer- 
ical difficulty of handling unbounded discontinuous solutions. In Fig. 6, the computed 
stationary solutions for (1) and their level lines are shown in the test example of Section 
4 (compare with Fig. 5). 

Despite the fact that the real sand flow is completely separated in the left and the 
right subregions of O, at the numerical level the flow travels through the grid points and 
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standing layer , N=31 rolling layer , N=31 




Figure 6. Numerical stationary solutions u and v of system (1) in the test 
example and their level lines 

then it can cross the separation line. More precisely, the transport path for sand from 

a point X in the right hand side should be the segment XP in Fig. 4; on the contrary, 
the algorithm splits this flow along vertical and horizontal segments connecting nodes and 
then part of this sand reaches the segment PQ even far from P (and from there eventually 
the left-hand side of the table). That is why the simple use at the discrete level of the 
same decomposition strategy adopted to characterize the stationary solutions is not able to 
reduce this phenomenon. As a test we applied in fact on the same uniform grid the scheme 
(33)-(37) separately in the two subregions of CI, with suitable wall boundary conditions 
on the cut (the PQ segment). The results (see Fig. 7) show an evident improvement of 
the solutions only in the left (that is the regular) subregion. 

Better results can be expected by coupling decomposition with suitable grid strategies. 
Keeping the uniform grid, the use of semi-lagrangian type schemes along characteristics 
should give a better trace of the correct transport directions. On the other hand, a different 
idea could be to employ unstructured grids (and mesh reflnements near the singularity 
regions) in order to improve the accuracy. The discussion of these approaches will be 
the goal of a forthcoming paper. The main difficulty is, anyway, that a sharp domain 
decomposition requires the a priori knowledge of the ridge set, which is not in general 
an easy task. For example, if we slightly modify the table in Fig. 4 by simply opening a 
symmetric portion of the boundary on its northern side ({(,T,y) : y = 1, 0.5 < x < 1}), 
the situation becomes completely different: a curved internal ridge appears, and with the 
help of the normal directions to the singular boundary points it subdivides the table into 
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5 10 15 20 25 30 5 10 15 20 25 30 

X X 



Figure 7. Numerical stationary solutions u and v of system (1) in the test 
example and their level lines: solution by decomposition 



standing layer u^. M= 31 

; rolling layerv^, %% 




Figure 8. Numerical stationary solutions u and v (view from above) of 
system (1) in the modified test example 



four distinct flow regions (sec Fig. 8, where the v surface is now seen from above, showing, 
in white, the ridge set profile). 
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